
function [c,ceq,Gc,Gceq]=eq_con_jac_thetas(pars,y,D,M,X,Z,R,sgi_nodes,sgi_weights,ep_sd,Sigma)

c=[];

[N,T,K]=size(X);
KZ=size(Z,3);
Kth=size(R,2);

M=reshape(M,N*T,1);
R=repmat(R,T,1);

alpha=repmat(pars(1:N,1),T,1);
theta=pars(N+(1:2*Kth),1);
theta=reshape(theta,Kth,2);
beta=pars(N+2*Kth+(1:2*N),1);
v=pars(3*N+2*Kth+(1:N),1);
f=repmat(pars(4*N+2*Kth+(1:N),1),T,1);
% vs=pars(5*N+2*Kth+(1:N),1);
xi=pars(6*N+2*Kth+(1:K),1);
gamma=pars(6*N+2*Kth+K+(1:KZ),1);
kappa=pars((6*N+2*Kth+K+KZ+1):end,1);

ep_sd=repmat(ep_sd,T,1);

if nargout<=2
    
    ZZ=zeros(N,N);
    for kz=1:KZ
        ZZ=ZZ+Z(:,:,kz)*gamma(kz,1);
    end 
    Q=diag(v.*ep_sd(1:N,1))*exp(-ZZ)*diag(v.*ep_sd(1:N,1));
    P=kron(eye(2),tril(Q)+tril(Q,-1)');
    P=P\eye(2*N);
    
    XX=zeros(N,T);
    Pd=[sqrt(diag(P\eye(2*N))),zeros(2*N,T-1)];
    B=[beta,zeros(2*N,T-1)];
    for t=1:T-1
        XX(:,t)=permute(X(:,t,:),[1,3,2])*xi;
        nm=find(M(N*(t-1)+(1:N),1)==1);
        DD=[diag(1-D(nm,t)),diag(D(nm,t))];
        P([nm;N+nm],[nm;N+nm])=P([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
        Pd(:,t+1)=sqrt(diag(P\eye(2*N)));
        beta([nm;N+nm],1)=beta([nm;N+nm],1)+P([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y(nm,t)-DD*beta([nm;N+nm],1))));
        B(:,t+1)=beta;
    end
    XX(:,T)=permute(X(:,T,:),[1,3,2])*xi;
    XX=reshape(XX,N*T,1);
    PdA=reshape(Pd(1:N,:),N*T,1);
    Pd=reshape(Pd(N+(1:N),:),N*T,1);
    BA=reshape(B(1:N,:),N*T,1);
    B=reshape(B(N+(1:N),:),N*T,1);

    UU0=exp(alpha+(R*theta(:,1)).*(PdA*sgi_nodes(:,1)'+BA+ep_sd*sgi_nodes(:,2)'));
    UU0=M.*(UU0./(1+UU0));
    UU1=exp(alpha+(R*theta(:,2)).*(Pd*sgi_nodes(:,1)'+B+ep_sd*sgi_nodes(:,2)')-f-XX-kappa);
    UU1=M.*(UU1./(1+UU1));

    ceq=(UU1-UU0)*sgi_weights;
    % ceq=[U1,U0];
    
else
    
    ZZ=zeros(N,N);
    for kz=1:KZ
        ZZ=ZZ+Z(:,:,kz)*gamma(kz);
    end 
    Q=diag(v.*ep_sd(1:N,1))*exp(-ZZ)*diag(v.*ep_sd(1:N,1));
    P=kron(eye(2),tril(Q)+tril(Q,-1)');
    P=P\eye(2*N);
    P=repmat(P,1,1,T);

    XX=zeros(N,T);
    Pd=[sqrt(diag(P(:,:,1)\eye(2*N))),zeros(2*N,T-1)];
    B=[beta,zeros(2*N,T-1)];
    for t=1:T-1
        XX(:,t)=permute(X(:,t,:),[1,3,2])*xi;
        P(:,:,t+1)=P(:,:,t);
        nm=find(M(N*(t-1)+(1:N),1)==1);
        DD=[diag(1-D(nm,t)),diag(D(nm,t))];
        P([nm;N+nm],[nm;N+nm],t+1)=P([nm;N+nm],[nm;N+nm],t+1)+DD'*(Sigma(nm,nm)\DD);
        Pd(:,t+1)=sqrt(diag(P(:,:,t+1)\eye(2*N)));
        beta([nm;N+nm],1)=beta([nm;N+nm],1)+P([nm;N+nm],[nm;N+nm],t+1)\(DD'*(Sigma(nm,nm)\(y(nm,t)-DD*beta([nm;N+nm],1))));
        B(:,t+1)=beta;
    end
    XX(:,T)=permute(X(:,T,:),[1,3,2])*xi;
    XX=reshape(XX,N*T,1);
    PdA=reshape(Pd(1:N,:),N*T,1);
    Pd=reshape(Pd(N+(1:N),:),N*T,1);
    BA=reshape(B(1:N,:),N*T,1);
    BD=reshape(B(N+(1:N),:),N*T,1);

    yA=PdA*sgi_nodes(:,1)'+BA+ep_sd*sgi_nodes(:,2)';
    yD=Pd*sgi_nodes(:,1)'+BD+ep_sd*sgi_nodes(:,2)';
    
    UU0=exp(alpha+(R*theta(:,1)).*yA);
    UU0=M.*(UU0./(1+UU0));
    UU1=exp(alpha+(R*theta(:,2)).*yD-f-XX-kappa);
    UU1=M.*(UU1./(1+UU1));

    ceq=(UU1-UU0)*sgi_weights;

    UU0=UU0.*(1-UU0);
    UU1=UU1.*(1-UU1);

    Gc=[];
    
    Gceq=zeros(N*T,length(pars));
    % alpha
    Gceq(:,1:N)=((UU1-UU0)*sgi_weights).*repmat(eye(N),T,1);
    % theta
    for k=1:Kth
        Gceq(:,N+k)=-(UU0.*R(:,k).*yA)*sgi_weights;
        Gceq(:,N+Kth+k)=(UU1.*R(:,k).*yD)*sgi_weights;
    end
    % beta, v, and gamma
    U1=UU1(1:N,:);
    U0=UU0(1:N,:);
    AA=eye(2*N);
    A=AA;
    DB=zeros(2*N,N+1); DS=DB;
    for n=1:N
          % beta
        Gceq(1:N,N+2*Kth+n)=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),zeros(2*N,1),A(:,n))*sgi_weights; %#ok<*SPRIX>
        Gceq(1:N,N+2*Kth+N+n)=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),zeros(2*N,1),A(:,N+n))*sgi_weights;
          % v
        [DS(:,n+1),DB(:,n+1)]=DB_vg(P(:,:,1),zeros(2*N,1),zeros(2*N,1),(1:N)',AA,Z,ZZ,v,ep_sd(1:N,1),n);
        Gceq(1:N,3*N+2*Kth+n)=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),DS(:,n+1),DB(:,n+1))*sgi_weights;
    end
      % gamma
    [DS(:,1),DB(:,1)]=DB_vg(P(:,:,1),zeros(2*N,1),zeros(2*N,1),(1:N)',AA,Z,ZZ,v,ep_sd(1:N,1),0);
    Gceq(1:N,6*N+2*Kth+K+(1:KZ))=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),DS(:,1),DB(:,1))*sgi_weights;
    for t=1:T-1
          % beta, v, and gamma
        U1=UU1(t*N+(1:N),:);
        U0=UU0(t*N+(1:N),:);
        nm=find(M(N*(t-1)+(1:N),1)==1);
        DD=[diag(1-D(nm,t)),diag(D(nm,t))];
        AA=eye(2*N);
        AA([nm;N+nm],[nm;N+nm])=AA([nm;N+nm],[nm;N+nm])-P([nm;N+nm],[nm;N+nm],t+1)\(DD'*(Sigma(nm,nm)\DD));
        A=AA*A;
        for n=1:N
              % beta
            Gceq(t*N+(1:N),N+2*Kth+n)=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),zeros(2*N,1),A(:,n))*sgi_weights;
            Gceq(t*N+(1:N),N+2*Kth+N+n)=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),zeros(2*N,1),A(:,N+n))*sgi_weights;
              % v
            [DS(:,n+1),DB(:,n+1)]=DB_vg(P(:,:,t+1),B(:,t+1)-B(:,t),DB(:,n+1),nm,AA,Z,ZZ,v,ep_sd(1:N,1),n);
            Gceq(t*N+(1:N),3*N+2*Kth+n)=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),DS(:,n+1),DB(:,n+1))*sgi_weights;
        end
          % gamma
        [DS(:,1),DB(:,1)]=DB_vg(P(:,:,t+1),B(:,t+1)-B(:,t),DB(:,1),nm,AA,Z,ZZ,v,ep_sd(1:N,1),0);
        Gceq(t*N+(1:N),6*N+2*Kth+K+(1:KZ))=J_bvg_thetas(U1,U0,R(1:N,:),theta,sgi_nodes(:,1),DS(:,1),DB(:,1))*sgi_weights;
    end    
    % f
    A=-(UU1)*sgi_weights;
    Gceq(:,4*N+2*Kth+(1:N))=A.*repmat(eye(N),T,1);
    % kappa
    Gceq(:,(6*N+2*Kth+K+KZ+1):end)=diag(A);
    % xi
    for k=1:K
        A=-(UU1.*reshape(X(:,:,k),N*T,1))*sgi_weights;
        Gceq(:,6*N+2*Kth+k)=A;
    end
    
    Gceq=sparse(Gceq');
    
end

